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INTRODUCTION 

Neutron stars are perhaps the most promising class of gravitational wave (GW) 
sources, and searches for such GW signals is particularly suited to the charac- 
teristics of the GEO600 detector (see Schutz "Getting Ready for GEO600 Data" 
gr-qc99 10033). However, the instantaneous GW frequency of such a source will 
evolve due to both intrinsic spindown effects and Doppler modulations induced by 
the motion of the Earth. Thus because of the large parameter space of likely signals, 
directly implemented optimal matched filtering is not computationally feasible. 

In response to this problem, Schutz and Papa have developed an alternative 
strategy: the Hough-Hierarchical search algorithm (see Schutz and Papa "End- 
to-End Algorithm for Hierarchical Area Searches for Long-Duration GW Sources 
for GEO600" gr-qc9905018). In order to carry out a blind search over a range of 
intrinsic GW frequencies, the following three stages must be calculated for each 
point in the parameter space of sky positions and intrinsic spindown parameters: 

Stage I: Calculate demodulated Fourier transforms (DeFTs) on an intermediate 
time baseline (of order 1 day) by combining FFTs of short durations (ap- 
proximately 30 minutes) of the time series data. In this context demodulated 
means that if there is a source at the sky position in question, and with the 
intrinsic spindown parameters in question, then all spindown and modulatory 
effects will have been correctly removed from the DeFTs: all signal power 
will be confined to one and the same frequency bin in each DeFT. This fre- 
quency is the intrinsic frequency of the source measured at the start of the 
observing time. It is expected that the total observing time will be of order 4 
months, and thus roughly 120 of these DeFTs will be calculated for each point 
in parameter space. 

Stage II: In general source parameters will not coincide exactly with those 
searched for, and residual frequency evolution and modulation will remain 
in the DeFTs. Thus, the peak in power associated with a given source may 



change frequency bins from DeFT to DeFT. Because of the relatively small 
time baseline of these DeFTs and the resultant poor signal-to-noise of any 
expected continuous GW signal, this evolution will be not directly apparent 
in the DeFTs, but can be recovered statistically using the Hough Transform 
algorithm. 

Stage III: Calculate DeFTs for candidate sources with the full frequency resolution 
of the total observation time, by combining the intermediate baseline DeFTs 
produced in stage I. 

Thus, during stage II, regions of the parameter space in which it is statistically 
unlikely that there are GW sources are eliminated from the search. Thereby, in 
stage III, the most computationally expensive part of the algorithm, the long time 
baseline DeFTs are calculated over only a very small fraction of parameter space 
and over a very small range of frequencies. 

In this paper we outline the methods used in the first and third stages of this 
algorithm in constructing a longer time baseline DeFT from a number of shorter 
time baseline FFTs or DeFTs. 

THE METHOD 

Consider a time series x a of total duration T, which has been divided into M 
short time series, each having N data points. Then the DeFT for a signal with a 
time independent amplitude and phase 27r$ afe (A) is 
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where the time indices are related by Na + j = a, and b is a long time baseline 
frequency index. In the following discussion Latin indices j, k, I always sum over 
N, while Greek indices sum over M. Note that $ a ft(A) is dependent on a vector A 
of parameters which characterize the signal one is searching for. In searching for 
GW signals from neutron stars these will include intrinsic spindown parameters, 
and the position of the source in the sky. If x a k is the matrix formed by carrying 
out Fourier transforms along the short time index j in x a j , then equation 1 can be 
written as 
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where the product Q a (b, X)P a k(b, A) is defined by the terms in square brackets, and 
Q a (b, A) contains all parts of the square brackets independent of the short time 
index j and short frequency index k. 
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In equation 2 we have effectively re-written equation 1, a long time baseline 
DeFT in the time domain, as a sum (a index) of short time baseline DeFTs in the 
frequency domain (k index), where Q a (b, X)P a k(b, A) are these frequency domain 
filters. In the presence of stationary noise with a flat spectrum, equation 2 is the 
optimal detector. However, through applying various approximations, the detector 
can be made "acceptably sub-optimal", in the sense that only a small fraction of 
power from a signal is lost in comparison to the optimal case, while achieving vast 
savings in computational cost. 

To illustrate these mathematical approximations it is instructive to discuss equa- 
tion 2 for a specific case of $ aj fe(A): a linearly varying frequency model, i.e. in the 
continuum limit <3>(t, / , fo) = fot + fot 2 , where f and f are the intrinsic frequency 
and spindown of the source respectively, and t is time. In the case of an actual 
search for GW signals from pulsars, <& a jb(\) will not be so simple. However, this 
model is sufficiently complex to effectively demonstrate all of the approximations 
to be discussed here. 

In discrete form, f , f ) can be written as $^73/(7) = (/3+Ml)(Na+j)/NM+ 
^(Na + j) 2 /N 2 M 2 , where the long time baseline frequency index b = (3 + Ml. 
The chosen discretization of the spindown parameter f = 7/T 2 is not practically 
appropriate. However, in an actual search, a grid of points in spindown parameter 
space will be chosen to ensure an acceptable loss of power from unresolved signals. 
Thus, in the following discussion, only searches for resolved f parameters will be 
considered. 

Approximation 1: By Taylor expanding the model phase function $(t) about the 
middle of each short duration time series (i.e. about j = N/2) and discarding 
terms of order (j/N) 2 = t 2 and higher, in the limit N — > 00 the function 
P a k{fi,l,l) is Re P ak ((3,l,~j) = sine x and Im P ak ((5,l,-f) = (1 -cosx)/x. In 
the phase model considered here x = —2ix (f3/M + / + (2a + 1)7/M 2 — k) and 
Q a {(3, I, 7) = exp {-2m{a(3/M + a^/M 2 )}. 

Approximation 2: Consider the case where the short time baseline is chosen such 
that the instantaneous model frequency f(t) = <&(t, fo, fo, fo, ■ ■ ■) does not 
move by more than one short time baseline frequency bin over the duration of a 
short time baseline data set, i.e. in the model discussed here \ f \T/M < M/T. 
Then for a given a, the function P a kip, A) will be peaked in power about the 
model frequency averaged over the duration of time associated with the ath 
short data set, i.e. about x = (the first three terms in the above definition 
of x are the index of this average model frequency). Thus only a few terms 
around this model frequency will contribute significantly to the summation 
over k in equation 2. 

Approximation 3: The semi-periodic nature of P a k(b, A) means that this function 
can be efficiently evaluated from a look-up table of values containing the peri- 
odic parts, and three further operations: to calculate one instance of P a k{b, A) 
will require only 8 floating point operations. 



Approximation 4: If one approximates the model frequency parameter (3 in the 
calculation of P a k(fl,l,^) as a fixed value, for example with f3 = fio, equation 
2 can be calculated as an FFT, i.e. 
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where n term relates to approximation 2, P ak ((3,l,~f) is defined above, and 
for the phase model discussed here Q a (j) = exp {— 2ni (a 2r y/M 2 — 7 /AM 2 )}. 
Thus for values of (3 sufficiently near to fio, the loss in power due to this approx- 
imation will be small. To obtain £3/(7) for other values of /3, the calculation 
must be repeated using another /3 . 
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RESULTS AND DISCUSSION 

Numerical tests have shown that if one chooses 10% as an acceptable loss in 
power in comparison to the optimal case, then Nfft = 8 and nt er m = 16 are the 
preferred parameter combination, if the short time baseline T/M is chosen such 
that in the phase model discussed here \ f \T/M < M/T. If one decides that only 
a 5% loss in optimal power is acceptable, then this can be achieved with the same 
parameters, but choosing T/M such that \ fo\T/M < M/2T. 

The computational cost of calculating one DeFT in stage I in floating point 
operations is 

C DeFT - 5.3 X 10- ( -L-) (^Z) , (4) 

V300 Hz/ \1 day/ V 8 A 16 /' yJ 

where B is the bandwidth of the search. This is comparable to the computational 
cost of the corresponding steps in the Hierarchical Stack / Slide algorithm of Brady 
and Creighton ("Searching for Periodic Sources with LIGO: Hierarchical Searches" 
gr-qc98 12014). The Hough-Hierarchical search algorithm also has a number of 
computational advantages. To calculate a given bandwidth of a DeFT requires 
only the FFT data from this bandwidth and an additional small overlap. Thus the 
algorithm can be easily parallelized by distributing data and work by bandwidth; 
and no communication between processors is required. Also, the complete three 
stage algorithm can be arranged in such a way that once a bandwidth of FFT 
data is read from disk by a processor, all computation required on this data can be 
carried out while this data is held in memory, thus time spent reading data from 
disk is a negligible fraction of the total computational time: each processor will 
need to read roughly 40 Mb from disk once every two weeks. Furthermore, little 
additional memory is required as workspace for stages I and III: less than 100 kb. 

The GEO600 data analysis team are currently working on coding this algorithm 
in a computationally optimal manner, as well as integrating this with the Hough 
Transform part of the procedure. 



